%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% This program replicates Panel C of Table 2 in the paper 
%%% "Optimal Portfolio Choice with Fat Tails and Parameter
%%% Uncertainty", Journal of Financial and Quantitative Analysis
%%% (2024), Raymond Kan & Nathan Lassance 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Note: due to the use of a finite number of simulations, 
%%% the final results are random and will not correspond 
%%% exactly to those reported in the paper.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
close all

%%% Load data
load Dataset25SBTM.txt;
R = Dataset25SBTM;

%%% Define variables
nobs = 10^1;%We use nobs=10^4 to create Table 2 in the paper
nu = 8;
gam = 1;

%%% Exponential integral
Expfun = @(n,x) integral(@(t)exp(-x.*t)./(t.^n),1,inf);

%%% Population parameters
N = 25;
onevec = ones(N,1);
mu = mean(R)';
sigma = cov(R,1);
invsigma = inv(sigma);
theta2 = mu'*invsigma*mu;
theta2g = (onevec'*invsigma*mu)^2/(onevec'*invsigma*onevec);
psi2 = theta2 - theta2g;
mug = (onevec'*invsigma*mu)/(onevec'*invsigma*onevec);
varg = mug^2/theta2g;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Known combination coefficients
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nobsk = 10^4;
%%% T = 60
T = 60;
rho = N/T;
M = eye(T)-(ones(T,1)*ones(1,T))./T;
%Normal
cnT60 = ((T-N-1)*(T-N-4)/(T*(T-2)))*theta2/(theta2+N/T);
c1nT60 = ((T-N-1)*(T-N-4)/(T*(T-2)))*psi2/(psi2+N/T);
c2nT60 = ((T-N-1)*(T-N-4)/(T*(T-2)))*(N/T)/(psi2+N/T);
%Student t (asymptotic)
y = @(x) (nu-2).*rho.*x./(2*(1-rho));
fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
eta = fzero(fun,0.5*(1+nu/(nu-2)));
varphi = 2*eta^2*(1-rho)/(nu-eta*(nu-2));
ctasympT60 = (1-rho)^2*theta2/((varphi/eta)*theta2+rho);
c1tasympT60 = (1-rho)^2*psi2/((varphi/eta)*psi2+rho);
c2tasympT60 = (1-rho)^2*(eta/varphi)*rho/((varphi/eta)*psi2+rho);
%Student t (exact)
Expk1vec = zeros(nobsk,1);
Expk2vec = zeros(nobsk,1);
Expk3vec = zeros(nobsk,1);
for n = 1:nobsk
     lambda = diag(sqrt((nu-2)./chi2rnd(nu,[T,1])));
     B = lambda*M*lambda;
     oneslambda = lambda*ones(T,1);
     Y = randn(T,N);
     Mat1 = inv(Y'*B*Y);
     Mat2 = Mat1^2;
     Expk1vec(n,1) = trace(Mat1);
     Expk2vec(n,1) = trace(Mat2);
     Expk3vec(n,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
end
a1 = (T-N-2)/N;
a2 = a1*(T-N-1)*(T-N-4)/(T-2);
a3 = a2/T;
a4 = (T-N-1)*(T-N-4)/(T*(T-2));
k1 = a1*mean(Expk1vec);
k2 = a2*mean(Expk2vec);
k3 = a3*mean(Expk3vec);
ctexactT60 = a4*(k1*theta2)/(k2*theta2+k3*rho);
c1texactT60 = a4*(k1*psi2)/(k2*psi2+k3*rho);
c2texactT60 = a4*(k3/k2)*(k1*rho)/(k2*psi2+k3*rho);
EU2fnT60 = 12*(1/(2*gam))*(T/(T-N-2))*(2*cnT60*k1*theta2-...
           cnT60^2*(k2*theta2+k3*N/T)/a4);
EU2ftasympT60 = 12*(1/(2*gam))*(T/(T-N-2))*(2*ctasympT60*k1*theta2-...
           ctasympT60^2*(k2*theta2+k3*N/T)/a4);
EU2ftexactT60 = 12*(1/(2*gam))*(T/(T-N-2))*(2*ctexactT60*k1*theta2-...
           ctexactT60^2*(k2*theta2+k3*N/T)/a4);
EU3fnT60 = 12*(1/(2*gam))*(T/(T-N-2))*(2*c1nT60*k1*theta2+2*c2nT60*k1*theta2g-...
           c1nT60^2*(k2*theta2+k3*N/T)/a4-c2nT60^2*k2*theta2g/a4-...
           2*c1nT60*c2nT60*k2*theta2g/a4);
EU3ftasympT60 = 12*(1/(2*gam))*(T/(T-N-2))*(2*c1tasympT60*k1*theta2+2*c2tasympT60*...
           k1*theta2g-c1tasympT60^2*(k2*theta2+k3*N/T)/a4-c2tasympT60^2*k2*...
           theta2g/a4-2*c1tasympT60*c2tasympT60*k2*theta2g/a4);
EU3ftexactT60 = 12*(1/(2*gam))*(T/(T-N-2))*(2*c1texactT60*k1*theta2+2*c2texactT60*...
           k1*theta2g-c1texactT60^2*(k2*theta2+k3*N/T)/a4-c2texactT60^2*k2*...
           theta2g/a4-2*c1texactT60*c2texactT60*k2*theta2g/a4);

%%% T = 120
T = 120;
rho = N/T;
M = eye(T)-(ones(T,1)*ones(1,T))./T;
%Normal
cnT120 = ((T-N-1)*(T-N-4)/(T*(T-2)))*theta2/(theta2+N/T);
c1nT120 = ((T-N-1)*(T-N-4)/(T*(T-2)))*psi2/(psi2+N/T);
c2nT120 = ((T-N-1)*(T-N-4)/(T*(T-2)))*(N/T)/(psi2+N/T);
%Student t (asymptotic)
y = @(x) (nu-2).*rho.*x./(2*(1-rho));
fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
eta = fzero(fun,0.5*(1+nu/(nu-2)));
varphi = 2*eta^2*(1-rho)/(nu-eta*(nu-2));
ctasympT120 = (1-rho)^2*theta2/((varphi/eta)*theta2+rho);
c1tasympT120 = (1-rho)^2*psi2/((varphi/eta)*psi2+rho);
c2tasympT120 = (1-rho)^2*(eta/varphi)*rho/((varphi/eta)*psi2+rho);
%Student t (exact)
Expk1vec = zeros(nobsk,1);
Expk2vec = zeros(nobsk,1);
Expk3vec = zeros(nobsk,1);
for n = 1:nobsk
     lambda = diag(sqrt((nu-2)./chi2rnd(nu,[T,1])));
     B = lambda*M*lambda;
     oneslambda = lambda*ones(T,1);
     Y = randn(T,N);
     Mat1 = inv(Y'*B*Y);
     Mat2 = Mat1^2;
     Expk1vec(n,1) = trace(Mat1);
     Expk2vec(n,1) = trace(Mat2);
     Expk3vec(n,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
end
a1 = (T-N-2)/N;
a2 = a1*(T-N-1)*(T-N-4)/(T-2);
a3 = a2/T;
a4 = (T-N-1)*(T-N-4)/(T*(T-2));
k1 = a1*mean(Expk1vec);
k2 = a2*mean(Expk2vec);
k3 = a3*mean(Expk3vec);
ctexactT120 = a4*(k1*theta2)/(k2*theta2+k3*rho);
c1texactT120 = a4*(k1*psi2)/(k2*psi2+k3*rho);
c2texactT120 = a4*(k3/k2)*(k1*rho)/(k2*psi2+k3*rho);
EU2fnT120 = 12*(1/(2*gam))*(T/(T-N-2))*(2*cnT120*k1*theta2-...
           cnT120^2*(k2*theta2+k3*N/T)/a4);
EU2ftasympT120 = 12*(1/(2*gam))*(T/(T-N-2))*(2*ctasympT120*k1*theta2-...
           ctasympT120^2*(k2*theta2+k3*N/T)/a4);
EU2ftexactT120 = 12*(1/(2*gam))*(T/(T-N-2))*(2*ctexactT120*k1*theta2-...
           ctexactT120^2*(k2*theta2+k3*N/T)/a4);
EU3fnT120 = 12*(1/(2*gam))*(T/(T-N-2))*(2*c1nT120*k1*theta2+2*c2nT120*k1*theta2g-...
           c1nT120^2*(k2*theta2+k3*N/T)/a4-c2nT120^2*k2*theta2g/a4-...
           2*c1nT120*c2nT120*k2*theta2g/a4);
EU3ftasympT120 = 12*(1/(2*gam))*(T/(T-N-2))*(2*c1tasympT120*k1*theta2+2*c2tasympT120*...
           k1*theta2g-c1tasympT120^2*(k2*theta2+k3*N/T)/a4-c2tasympT120^2*k2*...
           theta2g/a4-2*c1tasympT120*c2tasympT120*k2*theta2g/a4);
EU3ftexactT120 = 12*(1/(2*gam))*(T/(T-N-2))*(2*c1texactT120*k1*theta2+2*c2texactT120*...
           k1*theta2g-c1texactT120^2*(k2*theta2+k3*N/T)/a4-c2texactT120^2*k2*...
           theta2g/a4-2*c1texactT120*c2texactT120*k2*theta2g/a4);

%%% T = 240
T = 240;
rho = N/T;
M = eye(T)-(ones(T,1)*ones(1,T))./T;
%Normal
cnT240 = ((T-N-1)*(T-N-4)/(T*(T-2)))*theta2/(theta2+N/T);
c1nT240 = ((T-N-1)*(T-N-4)/(T*(T-2)))*psi2/(psi2+N/T);
c2nT240 = ((T-N-1)*(T-N-4)/(T*(T-2)))*(N/T)/(psi2+N/T);
%Student t (asymptotic)
y = @(x) (nu-2).*rho.*x./(2*(1-rho));
fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
eta = fzero(fun,0.5*(1+nu/(nu-2)));
varphi = 2*eta^2*(1-rho)/(nu-eta*(nu-2));
ctasympT240 = (1-rho)^2*theta2/((varphi/eta)*theta2+rho);
c1tasympT240 = (1-rho)^2*psi2/((varphi/eta)*psi2+rho);
c2tasympT240 = (1-rho)^2*(eta/varphi)*rho/((varphi/eta)*psi2+rho);
%Student t (exact)
Expk1vec = zeros(nobsk,1);
Expk2vec = zeros(nobsk,1);
Expk3vec = zeros(nobsk,1);
for n = 1:nobsk
     lambda = diag(sqrt((nu-2)./chi2rnd(nu,[T,1])));
     B = lambda*M*lambda;
     oneslambda = lambda*ones(T,1);
     Y = randn(T,N);
     Mat1 = inv(Y'*B*Y);
     Mat2 = Mat1^2;
     Expk1vec(n,1) = trace(Mat1);
     Expk2vec(n,1) = trace(Mat2);
     Expk3vec(n,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
end
a1 = (T-N-2)/N;
a2 = a1*(T-N-1)*(T-N-4)/(T-2);
a3 = a2/T;
a4 = (T-N-1)*(T-N-4)/(T*(T-2));
k1 = a1*mean(Expk1vec);
k2 = a2*mean(Expk2vec);
k3 = a3*mean(Expk3vec);
ctexactT240 = a4*(k1*theta2)/(k2*theta2+k3*rho);
c1texactT240 = a4*(k1*psi2)/(k2*psi2+k3*rho);
c2texactT240 = a4*(k3/k2)*(k1*rho)/(k2*psi2+k3*rho);
EU2fnT240 = 12*(1/(2*gam))*(T/(T-N-2))*(2*cnT240*k1*theta2-...
           cnT240^2*(k2*theta2+k3*N/T)/a4);
EU2ftasympT240 = 12*(1/(2*gam))*(T/(T-N-2))*(2*ctasympT240*k1*theta2-...
           ctasympT240^2*(k2*theta2+k3*N/T)/a4);
EU2ftexactT240 = 12*(1/(2*gam))*(T/(T-N-2))*(2*ctexactT240*k1*theta2-...
           ctexactT240^2*(k2*theta2+k3*N/T)/a4);
EU3fnT240 = 12*(1/(2*gam))*(T/(T-N-2))*(2*c1nT240*k1*theta2+2*c2nT240*k1*theta2g-...
           c1nT240^2*(k2*theta2+k3*N/T)/a4-c2nT240^2*k2*theta2g/a4-...
           2*c1nT240*c2nT240*k2*theta2g/a4);
EU3ftasympT240 = 12*(1/(2*gam))*(T/(T-N-2))*(2*c1tasympT240*k1*theta2+2*c2tasympT240*...
           k1*theta2g-c1tasympT240^2*(k2*theta2+k3*N/T)/a4-c2tasympT240^2*k2*...
           theta2g/a4-2*c1tasympT240*c2tasympT240*k2*theta2g/a4);
EU3ftexactT240 = 12*(1/(2*gam))*(T/(T-N-2))*(2*c1texactT240*k1*theta2+2*c2texactT240*...
           k1*theta2g-c1texactT240^2*(k2*theta2+k3*N/T)/a4-c2texactT240^2*k2*...
           theta2g/a4-2*c1texactT240*c2texactT240*k2*theta2g/a4);

Table2PanelCKnown = [EU2fnT60,EU2ftasympT60,EU2ftexactT60;...
     cnT60,ctasympT60,ctexactT60;EU2fnT120,EU2ftasympT120,EU2ftexactT120;...
     cnT120,ctasympT120,ctexactT120;EU2fnT240,EU2ftasympT240,EU2ftexactT240;...
     cnT240,ctasympT240,ctexactT240;EU3fnT60,EU3ftasympT60,EU3ftexactT60;...
     c1nT60,c1tasympT60,c1texactT60;c2nT60,c2tasympT60,c2texactT60;...
     EU3fnT120,EU3ftasympT120,EU3ftexactT120;c1nT120,c1tasympT120,c1texactT120;...
     c2nT120,c2tasympT120,c2texactT120;EU3fnT240,EU3ftasympT240,EU3ftexactT240;...
     c1nT240,c1tasympT240,c1texactT240;c2nT240,c2tasympT240,c2texactT240]


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Known combination coefficients
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% T = 60
T = 60;
rho = N/T;
M = eye(T)-(ones(T,1)*ones(1,T))./T;
a1 = (T-N-2)/N;
a2 = a1*(T-N-1)*(T-N-4)/(T-2);
a3 = a2/T;
a4 = ((T-N-1)*(T-N-4)/(T*(T-2)));
%Initialize vectors (c for two-fund, (c1,c2) for three-fund) 
%Six calibrations: 1) Normal (exact), 2) t (asymp), 3) t (exact),
%4) Elliptical (asymp), 5) Elliptical (exact), 6) Cross-validation
cvec1T60 = zeros(nobs,1);
c1vec1T60 = zeros(nobs,1);
c2vec1T60 = zeros(nobs,1);
cvec2T60 = zeros(nobs,1);
c1vec2T60 = zeros(nobs,1);
c2vec2T60 = zeros(nobs,1);
cvec3T60 = zeros(nobs,1);
c1vec3T60 = zeros(nobs,1);
c2vec3T60 = zeros(nobs,1);
cvec4T60 = zeros(nobs,1);
c1vec4T60 = zeros(nobs,1);
c2vec4T60 = zeros(nobs,1);
cvec5T60 = zeros(nobs,1);
c1vec5T60 = zeros(nobs,1);
c2vec5T60 = zeros(nobs,1);
cvec6T60 = zeros(nobs,1);
c1vec6T60 = zeros(nobs,1);
c2vec6T60 = zeros(nobs,1);
EU2f1T60 = zeros(nobs,1);
EU3f1T60 = zeros(nobs,1);
EU2f2T60 = zeros(nobs,1);
EU3f2T60 = zeros(nobs,1);
EU2f3T60 = zeros(nobs,1);
EU3f3T60 = zeros(nobs,1);
EU2f4T60 = zeros(nobs,1);
EU3f4T60 = zeros(nobs,1);
EU2f5T60 = zeros(nobs,1);
EU3f5T60 = zeros(nobs,1);
EU2f6T60 = zeros(nobs,1);
EU3f6T60 = zeros(nobs,1);
%Monte Carlo  
for n = 1:nobs

    if mod(n,nobs/100) == 0
        SimulationsLeftT60 = nobs-n
    end

    %Simulate returns
    lambda = diag(sqrt((nu-2)./chi2rnd(nu,1,T)));
    Y = randn(T,N);
    X = (mu + sqrtm(sigma)*Y'*lambda)';

    %Inputs
    muhat = mean(X)';
    sigmahat = cov(X,1);
    invsigmahat = inv(sigmahat);
    theta2hat = muhat'*invsigmahat*muhat;
    theta2hat = ((T-N-2)*theta2hat-N)/T + ((2*(theta2hat)^(N/2)...
                *((1+theta2hat)^(-(T-2)/2)))/(T*beta(N/2,(T-N)/2)*betainc(theta2hat/(1+theta2hat),N/2,(T-N)/2)));
    wghat = invsigmahat*onevec./(onevec'*invsigmahat*onevec);
    mughat = wghat'*muhat;
    varghat = wghat'*sigmahat*wghat;
    psi2hat = muhat'*invsigmahat*muhat - mughat^2/varghat;
    psi2hat = ((T-N-1)*psi2hat-(N-1))/T + ((2*(psi2hat)^((N-1)/2)...
              *((1+psi2hat)^(-(T-2)/2)))/(T*beta((N-1)/2,(T-N+1)/2)*betainc(psi2hat/(1+psi2hat),(N-1)/2,(T-N+1)/2)));
    [~,~,nuhat] = mlstudent(X);
    nuhat = min(1000,nuhat);
    tauhat = sum((X-muhat').^2,2);
    tauhat = tauhat./mean(tauhat);

    %1) Normal (exact)
    cvec1T60(n,1) = a4*theta2hat/(theta2hat+N/T);
    w = (cvec1T60(n,1)/gam)*invsigmahat*muhat;
    EU2f1T60(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
    c1vec1T60(n,1) = a4*psi2hat/(psi2hat+N/T);
    c2vec1T60(n,1) = a4*(N/T)/(psi2hat+N/T);
    w = (c1vec1T60(n,1)/gam)*invsigmahat*muhat + (c2vec1T60(n,1)/gam)*mughat*invsigmahat*onevec;
    EU3f1T60(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);

    %2) t (asymptotic)
    y = @(x) (nuhat-2).*rho.*x./(2*(1-rho));
    fun = @(x) y(x).*exp(y(x)).*Expfun(nuhat/2,y(x))-rho;
    eta = fzero(fun,0.5*(1+nuhat/(nuhat-2)));
    varphi = 2*eta^2*(1-rho)/(nuhat-eta*(nuhat-2));
    cvec2T60(n,1) = (1-rho)^2*theta2hat/((varphi/eta)*theta2hat+rho);
    w = (cvec2T60(n,1)/gam)*invsigmahat*muhat;
    EU2f2T60(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
    c1vec2T60(n,1) = (1-rho)^2*psi2hat/((varphi/eta)*psi2hat+rho);
    c2vec2T60(n,1) = (1-rho)^2*(eta/varphi)*rho/((varphi/eta)*psi2hat+rho);
    w = (c1vec2T60(n,1)/gam)*invsigmahat*muhat + (c2vec2T60(n,1)/gam)*mughat*invsigmahat*onevec;
    EU3f2T60(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);

    %3) t (exact)
    nobsk = 10^3;
    Expk1vec = zeros(nobsk,1);
    Expk2vec = zeros(nobsk,1);
    Expk3vec = zeros(nobsk,1);
    for m = 1:nobsk
         lambdahat = diag(sqrt((nuhat-2)./chi2rnd(nuhat,[T,1])));
         B = lambdahat*M*lambdahat;
         oneslambda = lambdahat*ones(T,1);
         Y = randn(T,N);
         Mat1 = inv(Y'*B*Y);
         Mat2 = Mat1^2;
         Expk1vec(m,1) = trace(Mat1);
         Expk2vec(m,1) = trace(Mat2);
         Expk3vec(m,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    k1 = a1*mean(Expk1vec);
    k2 = a2*mean(Expk2vec);
    k3 = a3*mean(Expk3vec);  
    cvec3T60(n,1) = a4*k1*theta2hat/(k2*theta2hat+k3*N/T);
    w = (cvec3T60(n,1)/gam)*invsigmahat*muhat;
    EU2f3T60(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
    c1vec3T60(n,1) = a4*k1*psi2hat/(k2*psi2hat+k3*rho);
    c2vec3T60(n,1) = a4*(k3/k2)*k1*rho/(k2*psi2hat+k3*rho);
    w = (c1vec3T60(n,1)/gam)*invsigmahat*muhat + (c2vec3T60(n,1)/gam)*mughat*invsigmahat*onevec;
    EU3f3T60(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);

    %4) Elliptical (asymptotic)
    fun = @(x) sum(1./(T-N+x.*tauhat.*N))-1;
    eta = fzero(fun,[1,10]);
    varphi = (1-rho)./(1/eta^2-sum((N.*tauhat.^2)./(T-N+N.*eta.*tauhat).^2));
    cvec4T60(n,1) = (1-rho)^2*theta2hat/((varphi/eta)*theta2hat+rho);
    w = (cvec4T60(n,1)/gam)*invsigmahat*muhat;
    EU2f4T60(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
    c1vec4T60(n,1) = (1-rho)^2*psi2hat/((varphi/eta)*psi2hat+rho);
    c2vec4T60(n,1) = (1-rho)^2*(eta/varphi)*rho/((varphi/eta)*psi2hat+rho);
    w = (c1vec4T60(n,1)/gam)*invsigmahat*muhat + (c2vec4T60(n,1)/gam)*mughat*invsigmahat*onevec;
    EU3f4T60(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
    
    %5. Elliptical (exact)
    nobsk = 10^3;
    lambdahat = diag(sqrt(tauhat));
    B = lambdahat*M*lambdahat;
    oneslambda = lambdahat*ones(T,1);
    Expk1vec = zeros(nobsk,1);
    Expk2vec = zeros(nobsk,1);
    Expk3vec = zeros(nobsk,1);
    for m = 1:nobsk
         Y = randn(T,N);
         Mat1 = inv(Y'*B*Y);
         Mat2 = Mat1^2;
         Expk1vec(m,1) = trace(Mat1);
         Expk2vec(m,1) = trace(Mat2);
         Expk3vec(m,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    k1 = a1*mean(Expk1vec);
    k2 = a2*mean(Expk2vec);
    k3 = a3*mean(Expk3vec);
    cvec5T60(n,1) = a4*k1*theta2hat/(k2*theta2hat+k3*N/T);
    w = (cvec5T60(n,1)/gam)*invsigmahat*muhat;
    EU2f5T60(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
    c1vec5T60(n,1) = a4*k1*psi2hat/(k2*psi2hat+k3*rho);
    c2vec5T60(n,1) = a4*(k3/k2)*k1*rho/(k2*psi2hat+k3*rho);
    w = (c1vec5T60(n,1)/gam)*invsigmahat*muhat + (c2vec5T60(n,1)/gam)*mughat*invsigmahat*onevec;
    EU3f5T60(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);

    %6) Cross-validation 
    kfold = 5;
    %Two-fund
    dc = 0.01;
    ckvec = 0:dc:1;
    Uk = zeros(length(ckvec),1);
    i = 1;
    for c = ckvec   
         Pk = zeros(T,1);
         for k = 1:kfold
              Xk = X;
              Xk(1+(k-1)*T/kfold:k*T/kfold,:) = [];
              Xoosk = X(1+(k-1)*T/kfold:k*T/kfold,:);
              muk = mean(Xk)';
              sigmak = cov(Xk,1);
              wk = (c/gam)*inv(sigmak)*muk;
              Pk(1+(k-1)*T/kfold:k*T/kfold,1) = Xoosk*wk;
          end
          Uk(i,1) = mean(Pk) - (gam/2)*var(Pk);
          i = i+1;
     end
     [~,index] = max(Uk);
     cvec6T60(n,1) = (index-1)*dc;
     w = (cvec6T60(n,1)/gam)*invsigmahat*muhat;  
     EU2f6T60(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
     %Three-fund
     dc = 0.025;
     c1kvec = 0:dc:1;
     c2kvec = 0:dc:1;
     Uk = zeros(length(c1kvec),length(c2kvec));
     i = 1;
     j = 1;
     for c1 = c1kvec   
          for c2 = c2kvec
               Pk = zeros(T,1);
                for k = 1:kfold
                     Xk = X;
                     Xk(1+(k-1)*T/kfold:k*T/kfold,:) = [];
                     Xoosk = X(1+(k-1)*T/kfold:k*T/kfold,:);
                     muk = mean(Xk)';
                     sigmak = cov(Xk,1);
                     wgk = inv(sigmak)*onevec./sum(inv(sigmak)*onevec);
                     mugk = wgk'*muk;
                     wk = (c1/gam)*inv(sigmak)*muk + (c2/gam)*mugk*inv(sigmak)*onevec;
                     Pk(1+(k-1)*T/kfold:k*T/kfold,1) = Xoosk*wk;
                end
                Uk(i,j) = mean(Pk) - (gam/2)*var(Pk);
                j = j+1;
          end
      j = 1;
      i = i+1;
      end
      [maxValue,~] = max(Uk(:));
      [index1,index2] = find(Uk == maxValue);
      c1vec6T60(n,1) = c1kvec(1) + (index1-1)*dc;
      c2vec6T60(n,1) = c2kvec(1) + (index2-1)*dc;
      w = (c1vec6T60(n,1)/gam)*invsigmahat*muhat + (c2vec6T60(n,1)/gam)*mughat*invsigmahat*onevec;
      EU3f6T60(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
end

%%% T = 120
T = 120;
rho = N/T;
M = eye(T)-(ones(T,1)*ones(1,T))./T;
a1 = (T-N-2)/N;
a2 = a1*(T-N-1)*(T-N-4)/(T-2);
a3 = a2/T;
a4 = ((T-N-1)*(T-N-4)/(T*(T-2)));
%Initialize vectors (c for two-fund, (c1,c2) for three-fund) 
%Six calibrations: 1) Normal (exact), 2) t (asymp), 3) t (exact),
%4) Elliptical (asymp), 5) Elliptical (exact), 6) Cross-validation
cvec1T120 = zeros(nobs,1);
c1vec1T120 = zeros(nobs,1);
c2vec1T120 = zeros(nobs,1);
cvec2T120 = zeros(nobs,1);
c1vec2T120 = zeros(nobs,1);
c2vec2T120 = zeros(nobs,1);
cvec3T120 = zeros(nobs,1);
c1vec3T120 = zeros(nobs,1);
c2vec3T120 = zeros(nobs,1);
cvec4T120 = zeros(nobs,1);
c1vec4T120 = zeros(nobs,1);
c2vec4T120 = zeros(nobs,1);
cvec5T120 = zeros(nobs,1);
c1vec5T120 = zeros(nobs,1);
c2vec5T120 = zeros(nobs,1);
cvec6T120 = zeros(nobs,1);
c1vec6T120 = zeros(nobs,1);
c2vec6T120 = zeros(nobs,1);
EU2f1T120 = zeros(nobs,1);
EU3f1T120 = zeros(nobs,1);
EU2f2T120 = zeros(nobs,1);
EU3f2T120 = zeros(nobs,1);
EU2f3T120 = zeros(nobs,1);
EU3f3T120 = zeros(nobs,1);
EU2f4T120 = zeros(nobs,1);
EU3f4T120 = zeros(nobs,1);
EU2f5T120 = zeros(nobs,1);
EU3f5T120 = zeros(nobs,1);
EU2f6T120 = zeros(nobs,1);
EU3f6T120 = zeros(nobs,1);
%Monte Carlo  
for n = 1:nobs

    if mod(n,nobs/100) == 0
        SimulationsLeftT120 = nobs-n
    end

    %Simulate returns
    lambda = diag(sqrt((nu-2)./chi2rnd(nu,1,T)));
    Y = randn(T,N);
    X = (mu + sqrtm(sigma)*Y'*lambda)';

    %Inputs
    muhat = mean(X)';
    sigmahat = cov(X,1);
    invsigmahat = inv(sigmahat);
    theta2hat = muhat'*invsigmahat*muhat;
    theta2hat = ((T-N-2)*theta2hat-N)/T + ((2*(theta2hat)^(N/2)...
                *((1+theta2hat)^(-(T-2)/2)))/(T*beta(N/2,(T-N)/2)*betainc(theta2hat/(1+theta2hat),N/2,(T-N)/2)));
    wghat = invsigmahat*onevec./(onevec'*invsigmahat*onevec);
    mughat = wghat'*muhat;
    varghat = wghat'*sigmahat*wghat;
    psi2hat = muhat'*invsigmahat*muhat - mughat^2/varghat;
    psi2hat = ((T-N-1)*psi2hat-(N-1))/T + ((2*(psi2hat)^((N-1)/2)...
              *((1+psi2hat)^(-(T-2)/2)))/(T*beta((N-1)/2,(T-N+1)/2)*betainc(psi2hat/(1+psi2hat),(N-1)/2,(T-N+1)/2)));
    [~,~,nuhat] = mlstudent(X);
    nuhat = min(1000,nuhat);
    tauhat = sum((X-muhat').^2,2);
    tauhat = tauhat./mean(tauhat);

    %1) Normal (exact)
    cvec1T120(n,1) = a4*theta2hat/(theta2hat+N/T);
    w = (cvec1T120(n,1)/gam)*invsigmahat*muhat;
    EU2f1T120(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
    c1vec1T120(n,1) = a4*psi2hat/(psi2hat+N/T);
    c2vec1T120(n,1) = a4*(N/T)/(psi2hat+N/T);
    w = (c1vec1T120(n,1)/gam)*invsigmahat*muhat + (c2vec1T120(n,1)/gam)*mughat*invsigmahat*onevec;
    EU3f1T120(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);

    %2) t (asymptotic)
    y = @(x) (nuhat-2).*rho.*x./(2*(1-rho));
    fun = @(x) y(x).*exp(y(x)).*Expfun(nuhat/2,y(x))-rho;
    eta = fzero(fun,0.5*(1+nuhat/(nuhat-2)));
    varphi = 2*eta^2*(1-rho)/(nuhat-eta*(nuhat-2));
    cvec2T120(n,1) = (1-rho)^2*theta2hat/((varphi/eta)*theta2hat+rho);
    w = (cvec2T120(n,1)/gam)*invsigmahat*muhat;
    EU2f2T120(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
    c1vec2T120(n,1) = (1-rho)^2*psi2hat/((varphi/eta)*psi2hat+rho);
    c2vec2T120(n,1) = (1-rho)^2*(eta/varphi)*rho/((varphi/eta)*psi2hat+rho);
    w = (c1vec2T120(n,1)/gam)*invsigmahat*muhat + (c2vec2T120(n,1)/gam)*mughat*invsigmahat*onevec;
    EU3f2T120(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);

    %3) t (exact)
    nobsk = 10^3;
    Expk1vec = zeros(nobsk,1);
    Expk2vec = zeros(nobsk,1);
    Expk3vec = zeros(nobsk,1);
    for m = 1:nobsk
         lambdahat = diag(sqrt((nuhat-2)./chi2rnd(nuhat,[T,1])));
         B = lambdahat*M*lambdahat;
         oneslambda = lambdahat*ones(T,1);
         Y = randn(T,N);
         Mat1 = inv(Y'*B*Y);
         Mat2 = Mat1^2;
         Expk1vec(m,1) = trace(Mat1);
         Expk2vec(m,1) = trace(Mat2);
         Expk3vec(m,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    k1 = a1*mean(Expk1vec);
    k2 = a2*mean(Expk2vec);
    k3 = a3*mean(Expk3vec);  
    cvec3T120(n,1) = a4*k1*theta2hat/(k2*theta2hat+k3*N/T);
    w = (cvec3T120(n,1)/gam)*invsigmahat*muhat;
    EU2f3T120(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
    c1vec3T120(n,1) = a4*k1*psi2hat/(k2*psi2hat+k3*rho);
    c2vec3T120(n,1) = a4*(k3/k2)*k1*rho/(k2*psi2hat+k3*rho);
    w = (c1vec3T120(n,1)/gam)*invsigmahat*muhat + (c2vec3T120(n,1)/gam)*mughat*invsigmahat*onevec;
    EU3f3T120(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);

    %4) Elliptical (asymptotic)
    fun = @(x) sum(1./(T-N+x.*tauhat.*N))-1;
    eta = fzero(fun,[1,10]);
    varphi = (1-rho)./(1/eta^2-sum((N.*tauhat.^2)./(T-N+N.*eta.*tauhat).^2));
    cvec4T120(n,1) = (1-rho)^2*theta2hat/((varphi/eta)*theta2hat+rho);
    w = (cvec4T120(n,1)/gam)*invsigmahat*muhat;
    EU2f4T120(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
    c1vec4T120(n,1) = (1-rho)^2*psi2hat/((varphi/eta)*psi2hat+rho);
    c2vec4T120(n,1) = (1-rho)^2*(eta/varphi)*rho/((varphi/eta)*psi2hat+rho);
    w = (c1vec4T120(n,1)/gam)*invsigmahat*muhat + (c2vec4T120(n,1)/gam)*mughat*invsigmahat*onevec;
    EU3f4T120(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
    
    %%% 5. Elliptical (exact)
    nobsk = 10^3;
    lambdahat = diag(sqrt(tauhat));
    B = lambdahat*M*lambdahat;
    oneslambda = lambdahat*ones(T,1);
    Expk1vec = zeros(nobsk,1);
    Expk2vec = zeros(nobsk,1);
    Expk3vec = zeros(nobsk,1);
    for m = 1:nobsk
         Y = randn(T,N);
         Mat1 = inv(Y'*B*Y);
         Mat2 = Mat1^2;
         Expk1vec(m,1) = trace(Mat1);
         Expk2vec(m,1) = trace(Mat2);
         Expk3vec(m,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    k1 = a1*mean(Expk1vec);
    k2 = a2*mean(Expk2vec);
    k3 = a3*mean(Expk3vec);
    cvec5T120(n,1) = a4*k1*theta2hat/(k2*theta2hat+k3*N/T);
    w = (cvec5T120(n,1)/gam)*invsigmahat*muhat;
    EU2f5T120(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
    c1vec5T120(n,1) = a4*k1*psi2hat/(k2*psi2hat+k3*rho);
    c2vec5T120(n,1) = a4*(k3/k2)*k1*rho/(k2*psi2hat+k3*rho);
    w = (c1vec5T120(n,1)/gam)*invsigmahat*muhat + (c2vec5T120(n,1)/gam)*mughat*invsigmahat*onevec;
    EU3f5T120(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);

    %%% 6) Cross-validation 
    kfold = 5;
    %Two-fund
    dc = 0.01;
    ckvec = 0:dc:1;
    Uk = zeros(length(ckvec),1);
    i = 1;
    for c = ckvec   
         Pk = zeros(T,1);
         for k = 1:kfold
              Xk = X;
              Xk(1+(k-1)*T/kfold:k*T/kfold,:) = [];
              Xoosk = X(1+(k-1)*T/kfold:k*T/kfold,:);
              muk = mean(Xk)';
              sigmak = cov(Xk,1);
              wk = (c/gam)*inv(sigmak)*muk;
              Pk(1+(k-1)*T/kfold:k*T/kfold,1) = Xoosk*wk;
          end
          Uk(i,1) = mean(Pk) - (gam/2)*var(Pk);
          i = i+1;
     end
     [~,index] = max(Uk);
     cvec6T120(n,1) = (index-1)*dc;
     w = (cvec6T120(n,1)/gam)*invsigmahat*muhat;  
     EU2f6T120(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
     %Three-fund
     dc = 0.025;
     c1kvec = 0:dc:1;
     c2kvec = 0:dc:1;
     Uk = zeros(length(c1kvec),length(c2kvec));
     i = 1;
     j = 1;
     for c1 = c1kvec   
          for c2 = c2kvec
               Pk = zeros(T,1);
                for k = 1:kfold
                     Xk = X;
                     Xk(1+(k-1)*T/kfold:k*T/kfold,:) = [];
                     Xoosk = X(1+(k-1)*T/kfold:k*T/kfold,:);
                     muk = mean(Xk)';
                     sigmak = cov(Xk,1);
                     wgk = inv(sigmak)*onevec./sum(inv(sigmak)*onevec);
                     mugk = wgk'*muk;
                     wk = (c1/gam)*inv(sigmak)*muk + (c2/gam)*mugk*inv(sigmak)*onevec;
                     Pk(1+(k-1)*T/kfold:k*T/kfold,1) = Xoosk*wk;
                end
                Uk(i,j) = mean(Pk) - (gam/2)*var(Pk);
                j = j+1;
          end
      j = 1;
      i = i+1;
      end
      [maxValue,~] = max(Uk(:));
      [index1,index2] = find(Uk == maxValue);
      c1vec6T120(n,1) = c1kvec(1) + (index1-1)*dc;
      c2vec6T120(n,1) = c2kvec(1) + (index2-1)*dc;
      w = (c1vec6T120(n,1)/gam)*invsigmahat*muhat + (c2vec6T120(n,1)/gam)*mughat*invsigmahat*onevec;
      EU3f6T120(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
end

%%% T = 240
T = 240;
rho = N/T;
M = eye(T)-(ones(T,1)*ones(1,T))./T;
a1 = (T-N-2)/N;
a2 = a1*(T-N-1)*(T-N-4)/(T-2);
a3 = a2/T;
a4 = ((T-N-1)*(T-N-4)/(T*(T-2)));
%Initialize vectors (c for two-fund, (c1,c2) for three-fund) 
%Six calibrations: 1) Normal (exact), 2) t (asymp), 3) t (exact),
%4) Elliptical (asymp), 5) Elliptical (exact), 6) Cross-validation
cvec1T240 = zeros(nobs,1);
c1vec1T240 = zeros(nobs,1);
c2vec1T240 = zeros(nobs,1);
cvec2T240 = zeros(nobs,1);
c1vec2T240 = zeros(nobs,1);
c2vec2T240 = zeros(nobs,1);
cvec3T240 = zeros(nobs,1);
c1vec3T240 = zeros(nobs,1);
c2vec3T240 = zeros(nobs,1);
cvec4T240 = zeros(nobs,1);
c1vec4T240 = zeros(nobs,1);
c2vec4T240 = zeros(nobs,1);
cvec5T240 = zeros(nobs,1);
c1vec5T240 = zeros(nobs,1);
c2vec5T240 = zeros(nobs,1);
cvec6T240 = zeros(nobs,1);
c1vec6T240 = zeros(nobs,1);
c2vec6T240 = zeros(nobs,1);
EU2f1T240 = zeros(nobs,1);
EU3f1T240 = zeros(nobs,1);
EU2f2T240 = zeros(nobs,1);
EU3f2T240 = zeros(nobs,1);
EU2f3T240 = zeros(nobs,1);
EU3f3T240 = zeros(nobs,1);
EU2f4T240 = zeros(nobs,1);
EU3f4T240 = zeros(nobs,1);
EU2f5T240 = zeros(nobs,1);
EU3f5T240 = zeros(nobs,1);
EU2f6T240 = zeros(nobs,1);
EU3f6T240 = zeros(nobs,1);
%Monte Carlo  
for n = 1:nobs

    if mod(n,nobs/100) == 0
        SimulationsLeftT240 = nobs-n
    end

    %Simulate returns
    lambda = diag(sqrt((nu-2)./chi2rnd(nu,1,T)));
    Y = randn(T,N);
    X = (mu + sqrtm(sigma)*Y'*lambda)';

    %Inputs
    muhat = mean(X)';
    sigmahat = cov(X,1);
    invsigmahat = inv(sigmahat);
    theta2hat = muhat'*invsigmahat*muhat;
    theta2hat = ((T-N-2)*theta2hat-N)/T + ((2*(theta2hat)^(N/2)...
                *((1+theta2hat)^(-(T-2)/2)))/(T*beta(N/2,(T-N)/2)*betainc(theta2hat/(1+theta2hat),N/2,(T-N)/2)));
    wghat = invsigmahat*onevec./(onevec'*invsigmahat*onevec);
    mughat = wghat'*muhat;
    varghat = wghat'*sigmahat*wghat;
    psi2hat = muhat'*invsigmahat*muhat - mughat^2/varghat;
    psi2hat = ((T-N-1)*psi2hat-(N-1))/T + ((2*(psi2hat)^((N-1)/2)...
              *((1+psi2hat)^(-(T-2)/2)))/(T*beta((N-1)/2,(T-N+1)/2)*betainc(psi2hat/(1+psi2hat),(N-1)/2,(T-N+1)/2)));
    [~,~,nuhat] = mlstudent(X);
    nuhat = min(1000,nuhat);
    tauhat = sum((X-muhat').^2,2);
    tauhat = tauhat./mean(tauhat);

    %1) Normal (exact)
    cvec1T240(n,1) = a4*theta2hat/(theta2hat+N/T);
    w = (cvec1T240(n,1)/gam)*invsigmahat*muhat;
    EU2f1T240(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
    c1vec1T240(n,1) = a4*psi2hat/(psi2hat+N/T);
    c2vec1T240(n,1) = a4*(N/T)/(psi2hat+N/T);
    w = (c1vec1T240(n,1)/gam)*invsigmahat*muhat + (c2vec1T240(n,1)/gam)*mughat*invsigmahat*onevec;
    EU3f1T240(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);

    %2) t (asymptotic)
    y = @(x) (nuhat-2).*rho.*x./(2*(1-rho));
    fun = @(x) y(x).*exp(y(x)).*Expfun(nuhat/2,y(x))-rho;
    eta = fzero(fun,0.5*(1+nuhat/(nuhat-2)));
    varphi = 2*eta^2*(1-rho)/(nuhat-eta*(nuhat-2));
    cvec2T240(n,1) = (1-rho)^2*theta2hat/((varphi/eta)*theta2hat+rho);
    w = (cvec2T240(n,1)/gam)*invsigmahat*muhat;
    EU2f2T240(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
    c1vec2T240(n,1) = (1-rho)^2*psi2hat/((varphi/eta)*psi2hat+rho);
    c2vec2T240(n,1) = (1-rho)^2*(eta/varphi)*rho/((varphi/eta)*psi2hat+rho);
    w = (c1vec2T240(n,1)/gam)*invsigmahat*muhat + (c2vec2T240(n,1)/gam)*mughat*invsigmahat*onevec;
    EU3f2T240(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);

    %3) t (exact)
    nobsk = 10^3;
    Expk1vec = zeros(nobsk,1);
    Expk2vec = zeros(nobsk,1);
    Expk3vec = zeros(nobsk,1);
    for m = 1:nobsk
         lambdahat = diag(sqrt((nuhat-2)./chi2rnd(nuhat,[T,1])));
         B = lambdahat*M*lambdahat;
         oneslambda = lambdahat*ones(T,1);
         Y = randn(T,N);
         Mat1 = inv(Y'*B*Y);
         Mat2 = Mat1^2;
         Expk1vec(m,1) = trace(Mat1);
         Expk2vec(m,1) = trace(Mat2);
         Expk3vec(m,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    k1 = a1*mean(Expk1vec);
    k2 = a2*mean(Expk2vec);
    k3 = a3*mean(Expk3vec);  
    cvec3T240(n,1) = a4*k1*theta2hat/(k2*theta2hat+k3*N/T);
    w = (cvec3T240(n,1)/gam)*invsigmahat*muhat;
    EU2f3T240(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
    c1vec3T240(n,1) = a4*k1*psi2hat/(k2*psi2hat+k3*rho);
    c2vec3T240(n,1) = a4*(k3/k2)*k1*rho/(k2*psi2hat+k3*rho);
    w = (c1vec3T240(n,1)/gam)*invsigmahat*muhat + (c2vec3T240(n,1)/gam)*mughat*invsigmahat*onevec;
    EU3f3T240(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);

    %4) Elliptical (asymptotic)
    fun = @(x) sum(1./(T-N+x.*tauhat.*N))-1;
    eta = fzero(fun,[1,10]);
    varphi = (1-rho)./(1/eta^2-sum((N.*tauhat.^2)./(T-N+N.*eta.*tauhat).^2));
    cvec4T240(n,1) = (1-rho)^2*theta2hat/((varphi/eta)*theta2hat+rho);
    w = (cvec4T240(n,1)/gam)*invsigmahat*muhat;
    EU2f4T240(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
    c1vec4T240(n,1) = (1-rho)^2*psi2hat/((varphi/eta)*psi2hat+rho);
    c2vec4T240(n,1) = (1-rho)^2*(eta/varphi)*rho/((varphi/eta)*psi2hat+rho);
    w = (c1vec4T240(n,1)/gam)*invsigmahat*muhat + (c2vec4T240(n,1)/gam)*mughat*invsigmahat*onevec;
    EU3f4T240(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
    
    %%% 5. Elliptical (exact)
    nobsk = 10^3;
    lambdahat = diag(sqrt(tauhat));
    B = lambdahat*M*lambdahat;
    oneslambda = lambdahat*ones(T,1);
    Expk1vec = zeros(nobsk,1);
    Expk2vec = zeros(nobsk,1);
    Expk3vec = zeros(nobsk,1);
    for m = 1:nobsk
         Y = randn(T,N);
         Mat1 = inv(Y'*B*Y);
         Mat2 = Mat1^2;
         Expk1vec(m,1) = trace(Mat1);
         Expk2vec(m,1) = trace(Mat2);
         Expk3vec(m,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    k1 = a1*mean(Expk1vec);
    k2 = a2*mean(Expk2vec);
    k3 = a3*mean(Expk3vec);
    cvec5T240(n,1) = a4*k1*theta2hat/(k2*theta2hat+k3*N/T);
    w = (cvec5T240(n,1)/gam)*invsigmahat*muhat;
    EU2f5T240(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
    c1vec5T240(n,1) = a4*k1*psi2hat/(k2*psi2hat+k3*rho);
    c2vec5T240(n,1) = a4*(k3/k2)*k1*rho/(k2*psi2hat+k3*rho);
    w = (c1vec5T240(n,1)/gam)*invsigmahat*muhat + (c2vec5T240(n,1)/gam)*mughat*invsigmahat*onevec;
    EU3f5T240(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);

    %%% 6) Cross-validation 
    kfold = 5;
    %Two-fund
    dc = 0.01;
    ckvec = 0:dc:1;
    Uk = zeros(length(ckvec),1);
    i = 1;
    for c = ckvec   
         Pk = zeros(T,1);
         for k = 1:kfold
              Xk = X;
              Xk(1+(k-1)*T/kfold:k*T/kfold,:) = [];
              Xoosk = X(1+(k-1)*T/kfold:k*T/kfold,:);
              muk = mean(Xk)';
              sigmak = cov(Xk,1);
              wk = (c/gam)*inv(sigmak)*muk;
              Pk(1+(k-1)*T/kfold:k*T/kfold,1) = Xoosk*wk;
          end
          Uk(i,1) = mean(Pk) - (gam/2)*var(Pk);
          i = i+1;
     end
     [~,index] = max(Uk);
     cvec6T240(n,1) = (index-1)*dc;
     w = (cvec6T240(n,1)/gam)*invsigmahat*muhat;  
     EU2f6T240(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
     %Three-fund
     dc = 0.025;
     c1kvec = 0:dc:1;
     c2kvec = 0:dc:1;
     Uk = zeros(length(c1kvec),length(c2kvec));
     i = 1;
     j = 1;
     for c1 = c1kvec   
          for c2 = c2kvec
               Pk = zeros(T,1);
                for k = 1:kfold
                     Xk = X;
                     Xk(1+(k-1)*T/kfold:k*T/kfold,:) = [];
                     Xoosk = X(1+(k-1)*T/kfold:k*T/kfold,:);
                     muk = mean(Xk)';
                     sigmak = cov(Xk,1);
                     wgk = inv(sigmak)*onevec./sum(inv(sigmak)*onevec);
                     mugk = wgk'*muk;
                     wk = (c1/gam)*inv(sigmak)*muk + (c2/gam)*mugk*inv(sigmak)*onevec;
                     Pk(1+(k-1)*T/kfold:k*T/kfold,1) = Xoosk*wk;
                end
                Uk(i,j) = mean(Pk) - (gam/2)*var(Pk);
                j = j+1;
          end
      j = 1;
      i = i+1;
      end
      [maxValue,~] = max(Uk(:));
      [index1,index2] = find(Uk == maxValue);
      c1vec6T240(n,1) = c1kvec(1) + (index1-1)*dc;
      c2vec6T240(n,1) = c2kvec(1) + (index2-1)*dc;
      w = (c1vec6T240(n,1)/gam)*invsigmahat*muhat + (c2vec6T240(n,1)/gam)*mughat*invsigmahat*onevec;
      EU3f6T240(n,1) = 12*(w'*mu - (gam/2)*w'*sigma*w);
end

Table2PanelCEstimated = ...
[mean(EU2f1T60),mean(EU2f2T60),mean(EU2f3T60),mean(EU2f4T60),mean(EU2f5T60),...
mean(EU2f6T60);mean(cvec1T60),mean(cvec2T60),mean(cvec3T60),mean(cvec4T60),...
mean(cvec5T60),mean(cvec6T60);...
mean(EU2f1T120),mean(EU2f2T120),mean(EU2f3T120),mean(EU2f4T120),mean(EU2f5T120),...
mean(EU2f6T120);mean(cvec1T120),mean(cvec2T120),mean(cvec3T120),mean(cvec4T120),...
mean(cvec5T120),mean(cvec6T120);...
mean(EU2f1T240),mean(EU2f2T240),mean(EU2f3T240),mean(EU2f4T240),mean(EU2f5T240),...
mean(EU2f6T240);mean(cvec1T240),mean(cvec2T240),mean(cvec3T240),mean(cvec4T240),...
mean(cvec5T240),mean(cvec6T240);...
mean(EU3f1T60),mean(EU3f2T60),mean(EU3f3T60),mean(EU3f4T60),mean(EU3f5T60),...
mean(EU3f6T60);mean(c1vec1T60),mean(c1vec2T60),mean(c1vec3T60),mean(c1vec4T60),...
mean(c1vec5T60),mean(c1vec6T60);mean(c2vec1T60),mean(c2vec2T60),mean(c2vec3T60),...
mean(c2vec4T60),mean(c2vec5T60),mean(c2vec6T60);...
mean(EU3f1T120),mean(EU3f2T120),mean(EU3f3T120),mean(EU3f4T120),mean(EU3f5T120),...
mean(EU3f6T120);mean(c1vec1T120),mean(c1vec2T120),mean(c1vec3T120),mean(c1vec4T120),...
mean(c1vec5T120),mean(c1vec6T120);mean(c2vec1T120),mean(c2vec2T120),mean(c2vec3T120),...
mean(c2vec4T120),mean(c2vec5T120),mean(c2vec6T120);...
mean(EU3f1T240),mean(EU3f2T240),mean(EU3f3T240),mean(EU3f4T240),mean(EU3f5T240),...
mean(EU3f6T240);mean(c1vec1T240),mean(c1vec2T240),mean(c1vec3T240),mean(c1vec4T240),...
mean(c1vec5T240),mean(c1vec6T240);mean(c2vec1T240),mean(c2vec2T240),mean(c2vec3T240),...
mean(c2vec4T240),mean(c2vec5T240),mean(c2vec6T240)]

TablePanelCComplete = [Table2PanelCKnown,Table2PanelCEstimated]

